using Distributions,Printf,Random,Statistics

const M0=10000000
const n_vec=[10,100,Inf]
const Y_mat=[-2.5 2.5;-0.05 0.05]

bounds  = zeros(size(Y_mat,1),2*length(n_vec)) 
results = Array{Float64,2}(undef,size(Y_mat,1),4*length(n_vec))
rng     = Xoshiro(2024)
mu0     = randn(rng,M0,2)
for (i,Y) in enumerate(eachrow(Y_mat))
    for (j,n) in enumerate(n_vec)
        global results
	global draws
        mu1 = ones(M0,1)*Y'+mu0/sqrt(n)
        mu  = mu1[mu1[:,1].<mu1[:,2],:]
        M   = size(mu,1)

        # prior 1: uniform prior
        theta = mu[:,1]+rand(rng,M).*(mu[:,2]-mu[:,1])
        results[i,4*(j-1)+1] = mean(theta.>0.0)	

	# prior 2: truncated normal prior
	theta = rand.(rng,Truncated.(Normal(),mu[:,1],mu[:,2]))
	results[i,4*(j-1)+2] = mean(theta.>0.0) 

	# prior 3: Extremely 'informative' prior i
        d = Float64.(rand(rng,Bernoulli(0.99),M))
        theta = d.*mu[:,1]+(ones(M,1)-d).*mu[:,2]
	results[i,4*(j-1)+3] = mean(theta.>0.0)

	# prior 4: Extremely 'informative' prior ii
        d = Float64.(rand(rng,Bernoulli(0.01),M))
        theta = d.*mu[:,1]+(ones(M,1)-d).*mu[:,2]
	results[i,4*j] = mean(theta.>0.0)

        if j<length(n_vec)
           bounds[i,2*j-1] = mean(mu[:,1].>0.0) 
           bounds[i,2*j]   = mean(mu[:,2].>0.0)
        else
           if Y[2]<0.0
                bounds[i,2*j-1:2*j] .= 0.0
            elseif Y[1]>0.0
                bounds[i,2*j-1:2*j] .= 1.0
            elseif (Y[1]<0.0)&(Y[2]>0.0)
                bounds[i,2*j-1] = 0.0 
                bounds[i,2*j]   = 1.0 
            end
        end            
    end
end 

io = open("watson_probability_bounds.tex","w")
@printf(io,"\\documentclass[10pt]{article}\n")
@printf(io,"\\usepackage{rotating}\n")
@printf(io,"\\pagestyle{empty}\n")
@printf(io,"\\begin{document}\n")
@printf(io,"\\tiny\n")
@printf(io,"\\begin{sidewaystable}\n")
@printf(io,"\\vspace*{1em}\n\n")
@printf(io,"\\begin{tabular}{ll|ll|llll|ll|llll|ll|llll}\\hline\n")
@printf(io," & & \\multicolumn{6}{|c|}{\$n=100\$} & \\multicolumn{6}{|c|}{\$n=500\$} & \\multicolumn{6}{|c}{\$n=\\infty\$} \\\\ \\cline{3-20} \n")
@printf(io," & & \\multicolumn{2}{|c|}{identified} & \\multicolumn{4}{|c|}{prior} & \\multicolumn{2}{c}{identified} & \\multicolumn{4}{|c|}{prior} & \\multicolumn{2}{c}{identified} & \\multicolumn{4}{|c}{prior}  \\\\ \\cline{5-8} \\cline{11-14} \\cline{17-20} \n")
@printf(io," \$Y_{1}\$ & \$Y_{2}\$ & \\multicolumn{2}{|c|}{interval} & 1 & 2 & 3 & 4 & \\multicolumn{2}{|c|}{interval} & 1 & 2 & 3 & 4 & \\multicolumn{2}{|c|}{interval} & 1 & 2 & 3 & 4 \\\\ \\hline \n")
for i=1:size(Y_mat,1)
    @printf(io,"%3.2f & %3.2f ",Y_mat[i,1],Y_mat[i,2])
    for j=1:length(n_vec)
        @printf(io," & %3.2f & %3.2f ",bounds[i,2*j-1],bounds[i,2*j])
        @printf(io," & %3.2f & %3.2f & %3.2f & %3.2f ",results[i,4*(j-1)+1],results[i,4*(j-1)+2],results[i,4*(j-1)+3],results[i,4*j])
    end
    @printf(io,"\\\\ \\hline \n") 
end
@printf(io,"\\end{tabular}\n\n")
@printf(io,"\\end{sidewaystable}\n")
@printf(io,"\\end{document}\n")
close(io)

